Topological phase transition in a generalized Kane-Mele-Hubbard model: A combined 
Quantum Monte Carlo and Green's function study 
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We study a generalized Kane-Mele-Hubbard model with third-neighbor hopping, an interacting 
two-dimensional model with a topological phase transition as a function of third-neighbor hopping, 
by means of the determinant projector Quantum Monte Carlo (QMC) method. This technique is 
essentially numerically exact on models without a fermion sign problem, such as the one we consider. 
We determine the interaction-dependence of the Z2 topological insulator/trivial insulator phase 
boundary by calculating the Z2 invariants directly from the single-particle Green's function. The 
interactions push the phase boundary to larger values of third-neighbor hopping, thus stabilizing the 
topological phase. The observation of boundary shifting entirely stems from quantum fluctuations. 
We also identify qualitative features of the single-particle Green's function which are computationally 
useful in numerical searches for topological phase transitions without the need to compute the full 
topological invariant. 
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PACS numbers: 71.10.Fd,71.70.Ej 

Introduction. -Ylecenibf, interest in a new state of mat- 
ter, topological insulators, has explodedii"— Z2 topo- 
logical insulators (TI) do not require interactions for 
their existence. However, intermediate strength electron- 
electron interactions have been shown to drive novel 
phases in slave-particle studies when the non-interacting 
limit is a TI Jr— Interaction effects have also appeared in 
experimental studies on the weakly correlated Bi-based 
rpjj^jM Moreover, a recently discovered Kondo topo- 
logical insulator— seems a promising venue to explore 
the strongly interacting limit. An essential challenge in 
many-body studies of TI systems is the direct character- 
ization of the interacting topological phases and phase 
transitions. This has largely been accomplished with ei- 
ther mean-field-like approaches or indirect evidence (such 
as the spontaneous appearance of an order parameter, or 
the closing of the single-particle excitation gap). Thus, 
it is important to understand the role of interactions 
in topological phases beyond the standard independent- 
particle and mean-field framework, ideally within an un- 
biased approach. 

Various approaches, including the entanglement 
entropy/spectruntiSiii or K- matrix theorjsi^ have also 
been proposed to characterize these topological phases. 
In the case of Z2 TI, topological invariants can 
be identified in terms of the single-particle Green's 
function i^^i^^'^° In certain cases, the frequency domain- 
winding-number— and a pole-expansion of the self- 
energySl could be useful in identifying interaction effects 
in a topological phase transition. The pole-structure of 
the Green's function in dynamical mean-field theory has 
been shown to be a powerful tool in the study of in- 
teraction effects in topological phases. The approach, 
however, still faces the limitation of being applicable only 
to local self-energy approximations. 



Interaction induced topological phase transitions 
have been studied in various models, including the 
Haldane-Hubbard model^^i^^ the Kane-Mele-Hubbard 
model— "—'^'^'^ and the interacting Bernevig-Hughes- 
Zhang model.—"— Within these models, there is also 
a topological phase transition at the single-particle 
level even without interaction)^^ which can be in- 
duced by a staggered onsite energy^ Rashba spin-orbit 
coupling,—!^ or a third-neighbor hopping, as we discuss 
in this Rapid Communication. To study this transition 
we use numerically exact determinant projector QMC to 
map out the interaction dependence of the topological 
phase transition as a function of third-neighbor hopping. 
We find that interactions tend to stabilize the topolog- 
ical phase, and we show the zero-frequency behavior of 
the Green's function as a function of third-neighbor hop- 
ping can be used to quantitatively determine the phase 
boundary. 

Model.-We consider the generalized Kane-Mele- 
Hubbard model (KMH) on the honeycomb lattice (unit 
cell sites labeled A and B) with real- valued third- 
neighbor hopping tsjv: H = Hq + Hjj with 

Ho = ^^'tCj^ + «^so Y^'^la^^J^J'y 



creates an electron 



and Hu = - 1)'. Here cj,, 

with spin-(T on site-i; the fermion number operator is 



Hi 



T ; a runs over 'f and ^. The spin-orbit cou- 



pling strength is Xso, and Vi 



-1 for counter-clockwise 



hopping with My = — 1 otherwise.— The spin-orbit cou- 
pHng term opens a bulk gap and drives the system to a 
Z2 TI for tsjv = 0,^ 



2 




FIG. 1: (Color online) (a) The first Brillouin zone of tlie lion- 
eycomb lattice. For tsN = 0, the Dirac points are located at 
Ki^2 = (± g^^ ,0) labeled by open and solid circles, respec- 
tively. The time-reversal invariant momentum (TRIM) points 
labeled by the green dots are F = (0,0), Mi, 2 = (±^, ^), 
and M3 = (0, 1^). (b) The noninteracting band structure 
of the generalized KM model Eq.(Il| at t^N = \t (here us- 
ing \so ~ 0.3t). Note that at the critical point separating 
the trivial insulator from the TI the Dirac cones shift to the 
TRIM points Mi ^2. 3, instead at Ki^2- 



The Brillouin zone (BZ) of the honeycomb lattice is 
shown in Fig. [l](a). For general t^N but vanishing A50, 
the model still exhibits a graphene-like band structure 
with gapless Dirac cones located at Ki,2 ~ (ib ^^ ,0), 
where a is the lattice constant. However, an arbitrary 
Xso will open a bulk gap and the generalized KM model 
turns into a Z2 TI or a trivial insulator depending on 
values of t^^- We find that for J7 = the critical value 
of the third- neighbor coupling tc is t^N = ^t. At tc, the 
bulk gap closes and the gapless Dirac cones shift away 
from the K points and move to time-reversal invariant 
momenta (TRIM), Afi,2 = j^) and M3 = (0, |f ). 

The band structure at the topological critical point is 
depicted in Fig. [1] (b). As t^N < ^t, the system is a 
Z2 TI, whereas as ^sat > if it is a trivial insulator. At 
the noninteracting level, the value of tc is independent of 
Xso- 

We next consider the Hubbard interaction Hu, given 
below Eq.([T]). In the presence of the Hubbard interaction, 
the topological phase boundary, tc-, shifts; a mean-field 
approach is unable to accurately determine tc for U 7^ 
0. In fact, we have verified that Hartree-Fock theory*'^ 
predicts no shift at all for U sufficiently small to avoid the 
magnetic transition. For U larger than this critical value, 
Uc, the topological band insulator state breaks down to a 
topologically trivial magnetic statej^"— 'Sli^i^'^ Since the 
generalized Kane-Mele Hubbard model we consider with 
the tsN term still preserves the essential band features of 
the Kane-Mele model, one can expect that in the strong 
coupling limit U > Uc our generalized model will also 
have a phase transition from the Z2 TI to the magnetic 
state. 

To study physics not captured within a mean-field the- 
ory, we choose a moderate Hubbard interaction U rela- 



tive to the bandwidth (small enough to avoid inducing 
the magnetic phase in the thermodynamic limit). Our 
main goal here is to demonstrate how the single-particle 
Green's functions computed within QMC in a fermion 
sign-free problem can be used to identify a correlated TI 
phase and topologically trivial insulating state. We leave 
a detailed analysis of the large U case for a future publica- 
tion. At half-filling, i.e., one fermion per site, the system 
has a particle-hole symmetry and the QMC simulations 
can perform accurate sampling without sign problems. 
Thus, one can accurately determine the phase boundary 
shifts at different U beyond the mean-field level. We find 
that as U increases, the critical value of t^N shifts to- 
wards a larger value, thus effectively stabilizing the Z2 
TI phase. 

Numerical Results.-ln our QMC calculations we use 
an imaginary time step At such that Art — 0.05 and 
an inverse temperature 8 such that Qt — 40. For the 
noninteracting case, for any finite Xso and at t^^ < tc, 
the system is a Z2 TI. We find that for Xso = 0-U, 
the model transitions to a magnetic state at U = 3t. 
To increase the threshold value of U needed to induce 
the magnetism, we consider a larger Xso = 0-4i or even 
Xso = t for different U. For comparison, in Fig. [2] we 
plot the value of the Z2 invariant as a function of t^j^ for 
different values of U. Open and solid symbols denote the 
noninteracting and interacting cases, respectively. Unless 
otherwise stated, we consider system sizes L x L = 6 x 6 
with periodic boundary conditions. We also study the 
finite size effects on the topological phase transition by 
comparing with 12 x 12 and 18 x 18 clusters. We find 
negligible changes in the transition point for these larger 
system sizes indicating that the location of the phase 
transition is already accurately captured in the L x L = 
6x6 system size. 

Using the single-particle Green's function we directly 
evaluate the Z2 invariant where 

(-ir= n (2) 

kiGTRIM 

and fj^. = {p^i\P\i2i) denotes the parity of the eigenstates 
of zero-frequency Green's functions^'^ (see details in sup- 
plementary information). Fig. [2] (a)-(c) depict the de- 
pendence of the Z2 invariant on t^N/t for U/t = 2, 3 
and 4. The open black squares denote the Z2 invariant 
given by tight-binding calculations with a 200 x 200 sys- 
tem size. The open red circles indicate the Z2 invariant 
calculated by QMC simulations for a 6 x 6 system at 
U = 0. The results are indistinguishable, confirming the 
accuracy of our QMC calculations in the non-interacting 
limit, and validating the 6x6 system size results. The 
location of topological phase boundary is tc — ^t. In 
the TI phase, only the Mi point is parity odd; the other 
three TRIM points are parity-even (i.e. ?7r = 3 — +1 
and fjMi = —1), so (—1)'' — —1. Across the transition 
upon increasing t^N, ffMi 2 3 change parity. F and Mi are 
parity-even whereas Af2,3 are parity-odd, so (—1)'' = 1. 
The blue solid triangles in Fig. [5] (a)-(c) depict the 



3 




0.25 0.30 0.35 0.40 0.45 



FIG. 2: (Color online) (a)-(c) Z2 invariant at U/t = 2, 3 
and 4 vs tsjv. The spin-orbital coupling is Xso ~ OAt. The 
black squares show the Z2 invariant given by the tight-binding 
calculations with 200 x 200. The red circle indicates the Z2 
invariant calculated by QMC simulations with 6 x 6 at 17 = 0. 
The blue solid triangles depict the Z2 invariant of the KMH 
model at [/ 7^ 0. (d)-(f) show the proportional coefficient 
determined by the relation: Go- (ki,0) = au,.a^ from QMC 
simulations vs Jsjv. All the open symbols indicate noninter- 
acting cases, i.e. (7 = 0. The solid symbols denote interacting 
cases. 



dependence of the Z2 invariant on is at for f7 ^ 0. With 
correlations, the parity properties of the TRIM points 
still remain and Eq. ([2]) to evaluate the Z2 invariant 
is still validj^i^ Strictly speaking, at each Monte Carlo 
measurement, the relation ryt = ±1 is not guaranteed. 
Hovifever, after a thousand QMC samplings, (?7k) = ±1 
with tiny numerical errors. At weak interaction, the 
phase boundary is barely seen to deviate. At U = 2t, the 
phase boundary is numerically estimated at t^N = 0.335t, 
which slightly deviates from tc = ^t. By increasing U, 



however, one can explicitly see that the interacting crit- 
ical points not only deviate from t/3 but move towards 
larger values, indicating the topological phase is stabi- 
lized by interactions. At U = 3t and 4t, the topolog- 
ical phase transitions take place at t^N — 0.341i and 
0.348i, respectively. Moreover, when Xso = t, the topo- 
logical phase boundary at [/ = 6i occurs at t^N = 0.352i. 
This indicates a significant (~ 10%) shift of the topolog- 
ical phase boundary driven by the Hubbard interaction. 
Moreover, no shift as a function of U is observed in a 
static Hartree-Fock mean-field approximation. It is thus 
the quantum fluctuations originating in the interactions 
that are important for shifting the phase boundary and 
stabilizing the topological phase. We believe this is likely 
to be a rather general result. 

Next, we investigate the single-particle Green's func- 
tion in our model. The parity operator is written 
as 1 and with inversion symmetry the Green's 

functions for each spin are simply proportional to a^: 
Gcr (ki,0) — aiciCr^ [or see Eq. (B2) in the supplemen- 
tal information]. In Fig. [D (d)-(f) we show the propor- 
tionality coefficient as a function of for finite U. 
For comparison, ak in the noninteracting case is also de- 
picted. At U = 0, we find the universal relations. 



a Ms and ^ ~aM2, 



(3) 



for all values of Xso and t^^. The values of ar behave 
smoothly as is varied through the topological critical 
points. However, the a coefficients on the other TRIM 
points are divergent at t^N — tc and change sign at a 
topological phase transition. At a critical point, the gap 
closes at the TRIM [c.f. Fig. [T] (b)] so the zero-frequency 
Green's functions are on the polesi^ Irrespective of the 
value of Xso, the location of the sign change is always at 
tc, consistent with the behavior of the Z2 invariant. 

Turning on the Hubbard interaction U, one can still 
observe the sign change in at at the topological phase 
transition. For finite U, the Green's functions retain 
their cr^-like form and the universal relations in Eq. ^ 
are still observed: — c^Aia ^^'^ — —0LM2 within 
QMC simulation errors, independent of the value of U jt. 
However, the positions of ak begin to change their signs 
away from t/3, as indicated by arrows in Figl2] (d)-(f), 
which label the topological phase boundaries in the in- 
teracting case. The locations for the sign change are 
consistent with the places where the Z^ invariants dra- 
matically jump. Note that at larger \] the magnitude of 
Q!k gradually vanish, but a sign change is still evident. 

Also in Figs. [2] (d)-(f) one can observe how the ak co- 
efhcients evolve upon increasing interactions. In the non- 
interacting case, the coefficients flip sign dramatically at 
tc = t/3. However, the values of ak decrease by increas- 
ing I] and the sign-flip behavior becomes more smooth 
with stronger interaction. This corresponds to a smeared 
phase boundary indicated by the Z^ invariant changes in 
FigsH] (a)-(c). Interestingly, away the topological phase 
transitions, e.g. tj,N = 0.2t and 0.5t, the coefficients ak 
for U ^ seem to return to their noninteracting values. 
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FIG. 3: (color online) The comparison of the single-particle 
Green's function coefficients a Mi 2 as a function of tsjv on 6 x 6 
(open symbols) and 12 x 12 (solid symbols) for (a) Aso = 0.4t 
and U = At (b) Ago = t and U = 6t. The insets indicate the 
comparison of the Z2 invariants vs t^N on the 6x6 and 12x12 
clusters with the same parameters. 
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FIG. 4: (Color online) Single-particle excitation gap Ac for 
different values of interaction: U /t = 0, 2, 3 and 4 with Ago = 
QAt on (a) 6 X 6 and (b) 12 x 12 clusters. For i7 / the 
single particle gap remains open across the topological phase 
transition, in contrast to the behavior in a non-interacting 
system. 



Therefore, interaction effects in ak are most apparent as 
Hn approaches the topological phase transition points. 

Finally, we investigate how finite size effects influence 
the topological phase transition boundaries with finite U . 
For this purpose, we compare the QMC results on 6 x 6 
and 12 x 12 in Fig. [31 For a comparison for generic pa- 
rameters, we consider at the Mi and M2 points for (a) 
\so = 0.4t and U and (b) \so = t and C/ = U. It 
is evident that while stronger interaction decreases aMi 
and in magnitude, the location of the sign change 
of ttk barely depends on the system size. Independent 
of system size, and switch sign at the same 

value ot t^N- Such behavior shows that the topological 
phase transition has a weak size dependence. The insets 
indicate the Z2 invariant for the two cases, also show- 
ing a small size dependence. However, on a small size, 
a stronger U [e.g. the inset of Fig. [3] (b)] will lead to 
a less sharp boundary determined by the Z2 invariants, 
compared to the ak behavior. For the same numerical 
accuracy, one can investigate the single-particle Green's 
functions on small sizes compared to the Z2 invariant 
to determine the topological phase transition boundary. 
This result implies that the single-particle Green's func- 
tion can be a powerful tool in detecting topological phase 
transitions in interacting systems without the need to 
evaluate the full topological invariant. (Although, this 
should certainly be checked in a few cases as it is the 
precise quantity that is used to distinguish the topologi- 
cal and non-topological phases.) 

We note that the single-particle excitation gap is not a 
reliable quantity to detect the topological phase bound- 
ary in finite-size interacting systems. The single-particle 
gap should close when undertaking the topological phase 
transition. As shown in Fig. 21 however, the single par- 
ticle gaps are finite at the phase transitions for t/ 7^ on 



the finite size simulations. Indeed, comparing the 6x6 
and 12 X 12 systems, we can clearly see the decay ten- 
dency upon increasing size. The QMC results on finite 
size scaling up to 18 x 18 confirm that around the phase 
boundaries the gaps vanish at L — >■ 00. Thus the behav- 
ior of the gaps is subject to strong finite size effect and 
the feature of vanishing excitation can be only observed 
in the thermodynamic limit. Moreover, the degree to 
which the phase transition is obscured by the single par- 
ticle gap also increases with increasing U . With 12 x 12, 
for U = 2t the gap seems to close around tc, but for 
U = At the behavior prevents one from determining the 
topological phase boundary. Therefore, in an interacting 
system, one should focus on the invariant itself and the 
single-particle Green's function. 

Summary. -'We have studied a generalized Kane-Mele- 
Hubbard model with an additional third-neighbor hop- 
ping term added. In the non-interacting limit the model 
exhibits a topological phase transition as a function of 
third-neighbor hopping. By choosing moderate Hubbard 
interactions without inducing antiferromagnetic order- 
ing, we study the topological phase transition in the 
interacting level. Using a numerically exact, fermion 
sign-free determinant projector QMC method, we have 
mapped the interaction-dependence of this phase bound- 
ary. Our main result is that interactions stabilize the 
topological phase by shifting the phase boundary to en- 
large the topological region. This effect is absent in a 
static Hartree-Fock mean-field theory, which indicates it 
is entirely the quantum fiuctuations associated with the 
interactions that enlarges the topological phase. We also 
show that the single-particle Green's function can more 
accurately determine the phase boundary than the Z2 in- 
variant (which is derived from it) for small system sizes. 
If this result can be reliably generalized, this could be a 
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useful insight in large-scale "numerical searches" for real 
materials with topological properties. The importance 
of fluctuation effects in our model also suggest that some 
density functional theory calculations could incorrectly 
predict the topological invariant of materials where quan- 
tum fluctuations are key to deciding the phase. 
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Appendix A: Sign-free determinant projector QMC 

The determinant QMC has been shown to be an excellent and unbiased approach to deal with strongly correlated 
system with Hubbard interactions 42riS In the projector algorithm, the ground state wave function |5'o) can be 
obtained using standard projection procedures on a trivial wave function |\E't), as long as one requires (^t|^o) 0. 
The expectation value of an observable A is obtained by 

(A) = lim ^ , (Al) 

The projection operator e~®^ can be discretized into many time slices e~^^ = [e^'^'^^]^^ with 9 = AtM where 
At <C 1 and M is the number of time slices with a large integer number; e^^^^ — e^'^'^(-^o+^'^' is the imaginary 
time-evolution propagator during At. The noninteracting ground state of is a good candidate for the trial wave 
function \^t)- With this trial wave function, we have confirmed that the determinant projector QMC is in a good 
agreement with our exact diagonalization results onaLxL = 3x3 system. By the first order Suzuki- Trotter 
decomposition, one can decompose e^'^'^^ as 

^-ArH ^^-ArHo^-ArHu^ (A2) 

where Hq is the single-particle Hamiltonian of the generalized Kane-Mele (KM) model as shown in Eq. (1) of the 
main text. Hu = ~ 1)^ involves 4 fermionic operators and cannot be represented in terms of single-particle 

basis. However, by the discrete S'i7(2)-invariant Hubbard-Stratonovich transformation^ the interacting imaginary 
time-evolution operator g-'^'^^u (^fgr JJ > 0) can be decomposed as 

^_Ar^(„,_i)2 ^ 1 j2 7(0e^VA^')(0(«.-i)+o(AT4), (A3) 

i=±l.±2 



where 7(±1) ^ I + V6/3, j{±2) = I - V6/3; r]{±l) = ±y 2(3 - VG) and 7y(±2) ±^^2(3 -I- ^6) are 4-component 
auxiliary fields determined by Monte Carlo samplings. Ref. (sTI - fssI provide pedagogical introductions about the QMC 
method. In this work, we employ Art = 0.05 in all the QMC simulations. 

In the determinant algorithm with the Suzuki- Trotter decomposition Eq. (|A2p and the Hubbard-Stratonovich 
transformation Eq. (IA3p . the denominator of Eq. (|AI[) reads as^^^^LM^^ (up to a constant factor) 

M M 

(*T|e"®^|*T) = {^T\Y[e-'^^"-\^T) ^ {'i'T\Y[e-^^"'>e-'^^""--\^T) (A4) 

r=l r=l 
M 

where "Ylii runs over possible auxiliary configurations rj^li^r), where i — 1^N,t = 1^ M; Hq is the matrix kernel 
of Hq with spin-(T. The probability weight p for an arbitrary auxiliary configuration {Tjili.r)} is simply denoted as^ 

p{M) = det (oMkr)]) det (oMkr)]) , (A5) 
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where det (o<,[?7(/,,^)]) = Trf^JJ^^^e'^''^-^^ e'V^v(h..)(n...-i)y when p < 0, QMC simulations 

meet notorious minus-sign problems. 

It has been proven that at half filling there exists a particle- hole symmetry in the KM model such that the probability 
is always positive-definitivej ^"^'^^ This character still remains even considering the real-valued third-neighbor hopping 
tsTV in the generalized KM model. The particle-hole transformation performs as 

where = —1 (^i = 1) if i belongs to A (B) sublattice. To show the positiveness of p{{i]}) in the generalized KM 
model, we employ the particle-hole transformation on the t^N hopping with but remain f unchanged. Upon such a 
transformation, the t^N tight-binding term turns out to be 

Note that the isjv hopping connects A and B lattices, so we have (— — 1. Therefore, upon the particle-hole 
transformation, Hq and Hq still have identical matrix kernels. 
The Hubbard interaction on J, transforms as 



^T—ri{h.r)[ni,i ~ -) 



which is the complex conjugate of Hu on f. Consequently, upon the particle-hole symmetry, one can have det(O^) = 
det(O^)* and the probability weight p = det(O^) det(O^) = |det(0-f)p being real positive. The QMC simulation in 
the half- filled generalized KM model is sign- free and numerically exact. 



Appendix B: single particle Green's functions 

Without sign problems, the QMC samplings provide highly accurate not only equal-time Green's functions but also 
time-displaced Green's functions in real space^i^ 

GAr,T) = {^o\cAr,r)cUO)\'^o), 

where r > 0. By performing double Fourier transformation we obtain the Green's functions in momentum space and 
Matsubar a frequency, i. e. Gcr(k, ia;„). 

It has been shown that zero frequency Green's functions are able to evaluate the Z2 invariant index in the interacting 
casei^ The Z2 invariant is determined by the parity of the eigenvectors of the inverse Green's functions 

[G(k„0)]-V.> = Ai.lM.>- 

Note that since there still exists an inversion symmetry in the generalized KMH model, the inverse Green's functions 
and the parity operator have simultaneous eigenvectors, i.e. f = Jy^jA'i)- In the (generalized) KM model, 
the parity operator exchanges A, B sublattices independent of spin index. Therefore, with the spinor convention 
^''1' = (c^^ ^ c^^), the parity operator is defined as P = / (g) cr"^ In the QMC simulations, the particle- 
hole symmetry provides G^(ki,0) — G^(ki,0) while kj is time-reversal invariant momentum (TRIM), i.e. k — — k. 
Therefore, we can directly diagonalize Go-(ki, 0) = [— i?k — S(ki, 0)]""'^ instead of inverse Green's functions for all k^ e 
TRIM points 

Gcr(ki,0)|/ii) = fiilfii), 

and choose the eigenvectors associated with positive eigenvalues {fii > 0, denoting occupied bands and are called 
right-zero^). In the honeycomb lattice, the TRIM points are F, Mi^2,3 as depicted in Fig. [5] Then we can employ 
the formalism proposed by Fu and Kane^^'^" to identify the Z2 invariant as 



(-ir- n ^^'^ 

kieTRIM 



(Bl) 
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FIG. 5: (Color online) The Brillouin zone of the honeycomb lattice. The time-reversal invariant momentum (TRIM) points 
labeled by the green dots are F — (0,0), Mi, 2 = (±-^^, j^), and M3 = (0, |^). The open and solid circles denote graphene 

Dirac points Ki,2 = (±^^,0) 



where 77^. — {ili\P\p,i). When v = for trivial insulator, whereas v = 1 is a, Z2 topological insulator. In the case of 
?7 = 0, fj^. = ±1. In the cases of finite J7, we find that {fi^.) — ±1 can be still obtained by sufficient QMC simulations. 
As t^N approaches the topological critical point, (—1)'^ will be smeared out and is laid between ±1. In this case, more 
QMC samplings are required for more accurate values. 

Note that since G^(ki,0) = G;(ki, 0), and G(ki,0) G-t-(ki, 0)® G;(ki, 0)] and P (= I®a'-^) have the simultaneous 
eigenvector sets, one has a relation: 



G<,(k,,0) = ak,cr^- 



(B2) 



In the context we show that in addition to the Z2 invariant, the proportional coefficient ak also plays another role to 
characterize a topological phase transition and even is more sensitive than v numerically. Upon the topological phase 
transition, the bulk gap will close at the TRIM points. Thus the single-particle Green's functions will be divergent 
on the polesi^ 

The relation Eq. (|B2|) should be expected both in the noninteracting and interacting cases. However, as ?7 ^ 
Eq. (|B2I) is not guaranteed in a single measurement in the QMC simulations. The proportionality relation between 
the zero-frequency Green's functions and the parity matrix can be recovered only upon enough samplings. To 
interpret this, we present the 6x6 benchmark results for the matrix elements of the zero-frequency Green's functions 
at k = Ml as a function of the number of measurements in Figs. [5] gij — [G(Mi, 0)]ij and m denotes the number 



(a)t3M=0-32t 



(b)t =0.37t 
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0.00 
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< llg,JI 
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1/Vm 
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).00 0.01 0.02 0.03 
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FIG. 6: (Color online) The matrix elements of the zero-frequency Green functions G(Mi,0) vs the number of samplings m at 
(a) tsjv = 0.32t and (b) t^M = Q.2)lt. Xso ~ 0-it and U = it. Ke[gij] and Im[(/,;j] denote the real part and imaginary part of 
[G(Mi,0)]ij, respectively; ||<?m|| denotes the diagonal component of G(Mi,0) in magnitudes. 



8 



CO 



0.10 



0.06 



0.02 



(a)^3„=0.4t, U=4t 

O t3„=0.34t 
A t =0.40t 



0.12 






(b) ^3o=0.4t, U=5t 


0.1 


■ t =0.30t 

3N 




O t3^=0.34t i . 


0.08 


-A t3^=0.40t 


0.06 




0.04 




0.02 








0.00 





0.0 0.1 0.2 0.3 0.4 
1/L 



0.20 



(/) 0.08 



0.04 



(c) >.3o=0.4t, U=6t 

O t3,=0.34t 
A t =0.40t 



S4 



0.0 0.1 0.2 0.3 0.4 

1/L 



0.0 0.1 0.2 0.3 0.4 

1/L 



0.12 


(d) X3„=t, U=6t 




0.10 








O t,„=0.35t 

3N 




0.08 


-A t =0 40t 




0.06 




- 


0.04 


• 




0.02 






0.00 







0.0 0.1 0.2 0.3 0.4 
1/L 



FIG. 7: (Color online) (a)-(c) The finite size scaling of the antiferromagnetic spin structure factor Saf/N vs 1/L at Xso ~ 0.4f 
and different U = 41, 5t, 6t. (d) Saf/N vs 1/L at Xso = t and U = 6t. Here N = 2x L^. 



of measurements. Ago = 0.4t and U ~ At are used. In this case the topological phase boundary is identified at 
tzN = 0.348t. We choose the value of ts^v close to the critical point. Fig. |6](a) shows t2,N = 0.32i in the Z2 topological 
insulator phase and (b) for tzN = 0.37/: in the trivial insulator. From the panels, it is evident that the structure of 
the Green's function does not fit Eq. (|B2[) as there are no sufficient samplings. At small to, the real parts of 1712 
and (721 arc not equal; furthermore 1712 and g2i have imaginary parts and gn and 322 are finite. However, one can 
see that upon sampling sufficient times Re[gi2] — Re[(72i], and meanwhile Im[(;i2], Im[(72i], ||5ii(22)|| go to zero. Thus 
in the to — > 00 limit, Eq. (jB2[) is recovered. Also note that q;mi = Re[(7i2] in both cases indicate opposite sign as 
observed by the signature of the topological phase transition. Moreover, by such to scaling, we also confirm that the 
value of the Z2 invariant also shows monotonically close to ±1. In our paper we choose the value of to large enough 
to determine the structure and extract the coefficients. 



Appendix C: Critical Hubbard interactions for antiferromagnetism 

In the generalized KMH model a strong Hubbard interaction can also derive the antiferromagnetic (AF) ordering, 
due to the bipartite lattice structure. Similarly to the KMH model (with t^M = 0) 54,62,63 ^ ^j-^g generalized KMH 
model, finite values of Ago also break the SU{2) symmetry down to the U{1) symmetry and the dominant magnetism 
behavior lies on x-y plane. The planar spin structure factor can be defined as^i^ 



(— 1)''' — f (— 1) for i € A{B) sublattice. This is similar to determining the Neel type ordering in a square lattice by 
using the antiferromagnetic spin structure factor at k = (tt, tt). To identify whether there exists the antiferromagnetism 
in the thermodynamic limit, we study the finite size scaling behavior of Saf at L — > 00. Generally speaking, the 
spin-orbital coupling will suppress AF ordering and larger A^o's are associated with larger UcS to induce the AF 
ordering. Note that due to the presence of third nearest neighboring hopping /37V which favors the Neel pattern in 
the second order perturbation, the threshold interaction Uc in the generalized KMH model is smaller than that in the 
KMH model. 

The QMC results on Saf/N vs 1/L are shown in Figs. [T] In (a), we can see that, for Xso ~ 0.4t, U — At is not 
sufficiently large to induce the AF ordering. At U — 5/, Saf is enhanced and the U value is close to the critical 
value to drive the AF ordering. In (c) under the interaction U = 6t, Saf saturates to finite values at 1/L 0, 
suggesting that the AF ordering exists in the thermodynamic limit. Fig. [7] (d) depicts the case of C/ = 6/ but at 
Xso = t- Compared to (c), where an AF ordering is induced, the structure factor in (d) still goes to zero in the 
L — >■ 00 limit. Thus, stronger spin-orbital couplings obviously suppress the existence of AF ordering and raise values 
of critical interactions in the generalized KMH model. 
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Appendix D: Single-particle excitation 

In this subsection, we present the approach to evaluate the single-particle excitation (charge gap) Ac in the QMC 
simulations. The charge gap is defined as the energy cost to add a particle into (or remove a particle from) the system 
composed of TV particles. Assuming that we have H\^^+'^) = E^^+^l^-^+i) and = then the charge 

gap reads Ac = Eq~^^ — Eq . It can be obtained via calculating the on-site time-displaced Green's functions which 
are written as 

G(f=0,T) = l^G.(z,i;r) 

Therefore, at large r, we have G{f = 0,r) ~ e"'^^" and then one can find the slope of In G'(r = 0,t) at large r 
to determine the value of Ac. Refs. PgUstMb^ and [s^ provide the detailed descriptions. The evaluation of the 
excitation by the on-site single-particle Green's function can determine the value of the single-particle excitation 
without concerning about specific momentum points, e.g. Ac(k). (Note that the gap of the KM model with A = 
closes at the Dirac points Ki 2, whereas the gap of the generalized KM model with closes at Mi 23.) 
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